###########################################################################
## Attributions of Fault and Support for Government Assistance for Workers
## Displaced by Technology
###########################################################################
## Seth Werfel
## Christopher Witko
## Tobias Heinrich
###########################################################################
## 3) Make graphs
###########################################################################



## Quick estimation of raw data, NP bootstrap
#############################################
load(file="output/Data_prepped.Rdata")
n_sim <- 1000
T_vec <- c("Domestic", "Foreign", "Tech")
prop <- array(NA, dim=c(n_sim, 2, 3, 7),
              dimnames=list(Iter=NULL, Outcome=c("Fault", "Support"),
                            Treat=T_vec, K=1:7))
for(i in 1:n_sim)
{
  for(m in 1:3)
  {
    this_data <- subset(data, Treat == T_vec[m])
    ind <- sample(1:nrow(this_data), size=nrow(this_data), replace=TRUE)

    for(j in 1:7)
    {
      prop[i, 1, m, j] <- weighted.mean(this_data[ind, "WFault"] == j, this_data[ind, "weight"])
      prop[i, 2, m, j] <- weighted.mean(this_data[ind, "Support"] == j, this_data[ind, "weight"])
    }
  }
}
dimnames(prop)[[3]][3] <- "Technology"
prop <- adply(.data=prop, .margins=2:4, 
              .fun=function(x) data.frame(Mean=mean(x), Q025=quantile(x, 0.025),
                                          Q975=quantile(x, 0.975)))


## Proportion of responses by treatment
g <- ggplot(data=subset(prop, Outcome == "Fault"), aes(x=K, ymin=Q025, ymax=Q975, y=Mean, group=factor(Treat), colour=factor(Treat)))
g <- g + theme_minimal()
g <- g + geom_hline(yintercept=0, size=.8)
g <- g + geom_point(size=3, position=position_dodge(width=.25))
g <- g + geom_linerange(size=0.8, position=position_dodge(width=.25))
g <- g + scale_x_discrete(breaks=1:7, labels=rev(c("Strongly\nagree", "Agree", "Somewhat\nagree", 
                                                   "Neither", "Somewhat\ndisagree", "Disagree", "Strongly\ndisagree")))
g <- g + xlab("Agree/ disagree that worker is at fault") + ylab("Proportion of responses")
g <- g + scale_colour_manual(values=c("black", "gray80", "gray50"), name="Treatment")
g <- g + theme(legend.position=c(0.9, 0.8))
g <- g + theme(axis.text=element_text(size=rel(1.0)))
ggsave(file="output/Proportion-Fault.pdf", width=11, height=4.5)


g <- ggplot(data=subset(prop, Outcome == "Support"), aes(x=K, ymin=Q025, ymax=Q975, y=Mean, group=factor(Treat), colour=factor(Treat)))
g <- g + theme_minimal()
g <- g + geom_hline(yintercept=0, size=.8)
g <- g + geom_point(size=3, position=position_dodge(width=.25))
g <- g + geom_linerange(size=0.8, position=position_dodge(width=.25))
g <- g + scale_x_discrete(breaks=1:7, labels=rev(c("Strongly\nsupport", "Support", "Somewhat\nsupport", 
                                                   "Neither", "Somewhat\noppose", "Oppose", "Strongly\noppose")))
g <- g + xlab("Support for worker getting unemployed benefits") + ylab("Proportion of responses")
g <- g + scale_colour_manual(values=c("black", "gray80", "gray50"), name="Treatment", guide="none")
g <- g + theme(axis.text=element_text(size=rel(1.0)))
ggsave(file="output/Proportion-Support.pdf", width=11, height=4.5)




## Load estimation results
##########################
load(file="output/Mediation_results.Rdata")
load(file="output/Main_difference.Rdata")
load(file="output/Main_predictions.Rdata")

## Some tweaks
store_med$Treat <- ifelse(store_med$Treat == "Foreign", "Foreign", "Technology")
store_med$Sig <- ifelse(store_med$PP >= 0.975, "Sig", ifelse(store_med$PP <= 0.025, "Sig", "NS"))
main_diff$Sig <- ifelse(main_diff$PP >= 0.975, "Sig", ifelse(main_diff$PP <= 0.025, "Sig", "NS"))
main_diff$Sig[main_diff$Mean == 0] <- "NS"



## Main treatment effects; levels
#################################
g <- ggplot(data=subset(main_diff, Y == "Fault"), aes(x=K, shape=Sig, ymin=Q025, ymax=Q975, y=Mean, group=factor(Treat), colour=factor(Treat)))
g <- g + theme_minimal()
g <- g + scale_shape_manual(values=c(1, 19), guide="none")
g <- g + geom_hline(yintercept=0, size=.8)
g <- g + geom_point(size=3, position=position_dodge(width=.25))
g <- g + geom_linerange(size=0.8, position=position_dodge(width=.25))
g <- g + scale_x_discrete(breaks=1:7, labels=rev(c("Strongly\nagree", "Agree", "Somewhat\nagree", 
                                                   "Neither", "Somewhat\ndisagree", "Disagree", "Strongly\ndisagree")))
g <- g + xlab("Agree/ disagree that worker is at fault") + ylab("Difference in probability")
g <- g + scale_colour_manual(values=c("black", "gray80", "gray50"), name="Treatment")
g <- g + theme(legend.position=c(0.9, 0.89))
g <- g + theme(strip.text=element_text(hjust=0.0, vjust=1, size=rel(0.82)),
               panel.spacing = unit(1, "lines"), 
               axis.line.x=element_line(size=.0000001),
               axis.text=element_text(size=rel(1.05)),
               plot.title=element_text(size=rel(1.3), face="bold"))
ggsave(file="output/Treatment-Fault.pdf", width=11, height=5)


g <- ggplot(data=subset(main_diff, Y == "Support"), aes(shape=Sig, x=K, ymin=Q025, ymax=Q975, y=Mean, group=factor(Treat), colour=factor(Treat)))
g <- g + scale_shape_manual(values=c(1, 19), guide="none")
g <- g + theme_minimal()
g <- g + geom_hline(yintercept=0, size=.8)
g <- g + geom_point(size=3, position=position_dodge(width=.25))
g <- g + geom_linerange(size=0.8, position=position_dodge(width=.25))
g <- g + scale_x_discrete(breaks=1:7, labels=rev(c("Strongly\nsupport", "Support", "Somewhat\nsupport", 
                                                   "Neither", "Somewhat\noppose", "Oppose", "Strongly\noppose")))
g <- g + xlab("Support for worker getting unemployed benefits") + ylab("Difference in probability")
g <- g + scale_colour_manual(values=c("black", "gray80", "gray50"), name="Treatment", guide="none")
g <- g + theme(axis.text=element_text(size=rel(1.05)))
ggsave(file="output/Treatment-Support.pdf", width=11, height=5)



## Main causal mediation graphs
################################
tmp <- subset(store_med, Comparison == "Main" & Estimand == "Indirect")
g <- ggplot(data=tmp, aes(shape=Sig, x=factor(Y), ymin=Q025, ymax=Q975, y=Mean))
g <- g + facet_wrap(~ Treat, ncol=2)
g <- g + scale_shape_manual(values=c(1, 19), guide="none")
g <- g + theme_minimal()
g <- g + geom_hline(yintercept=0, size=.8)
g <- g + geom_point(size=3, position=position_dodge(width=.25), colour="black")
g <- g + geom_linerange(size=0.8, position=position_dodge(width=.25), colour="black")
#g <- g + scale_colour_manual(values=c("black", "gray60"))
g <- g + theme(strip.text=element_text(hjust=0.0, vjust=1, face="bold", size=rel(1.15)),
               panel.spacing = unit(1, "lines"), 
               axis.line.x=element_line(size=.0000001),
               axis.text=element_text(size=rel(0.89)),
               plot.title=element_text(size=rel(1.3), face="bold"))
g <- g + scale_x_discrete(breaks=1:7, labels=rev(c("Strongly\nsupport", "Support", "Somewhat\nsupport", 
                                                   "Neither", "Somewhat\noppose", "Oppose", "Strongly\noppose")))
g <- g + xlab("Support for worker getting unemployed benefits") + ylab("Difference in probability")
g <- g + theme(legend.position=c(.9, .18))
ggsave(file="output/Mediation-Main.pdf", width=11, height=5)




## Heterogeneous comparison graphs
##################################
u_comparison <- unique(store_med$Comparison)
for(i in 1:length(u_comparison))
{
  tmp <- subset(store_med, Comparison == u_comparison[i])
  
  g <- ggplot(data=tmp, aes(x=Y, ymin=Q025, y=Mean, ymax=Q975, shape=Sig))
  g <- g + facet_grid(Estimand ~ Treat)
  g <- g + theme_minimal()
  g <- g + scale_shape_manual(values=c(1, 19), guide="none")
  g <- g + geom_hline(yintercept=0, size=.8)
  g <- g + geom_point(size=3, position=position_dodge(width=.25))
  g <- g + geom_linerange(size=0.8, position=position_dodge(width=.25))
  g <- g + scale_x_discrete(breaks=1:7, labels=rev(c("Strongly\nsupport", "Support", "Somewhat\nsupport", 
                                                     "Neither", "Somewhat\noppose", "Oppose", "Strongly\noppose")))
  g <- g + xlab("Support for worker getting unemployed benefits") + ylab("Difference in probability")
  g <- g + theme(strip.text=element_text(hjust=0.0, vjust=1, face="bold", size=rel(1.09)),
                 panel.spacing = unit(1, "lines"), 
                 axis.line.x=element_line(size=.0000001),
                 axis.text=element_text(size=rel(0.90)),
                 plot.title=element_text(size=rel(1.3), face="bold"))
  
  ggsave(filename=paste0("output/Heterogeneous_", tolower(tmp$Comparison[1]), ".pdf"), width=11, height=6)
}



## Garbage collection
#####################
rm(list=ls())



